Confiabilidade de Sistemas: Falhas e Reparo
30 de setembro de 2025
Do estado 2: - Permanecer em 2: (1-p)² - Ir para 1: 2p(1-p) - Ir para 0: p²
Do estado 1: - Ir para 2: (1-p)q - Permanecer em 1: pq + (1-p)(1-q) - Ir para 0: p(1-q)
Do estado 0: - Ir para 2: q² - Ir para 1: 2q(1-q) - Permanecer em 0: (1-q)²
\[ P = \begin{pmatrix} (1-q)^2 & 2q(1-q) & q²\\ p(1-q) & pq + (1-p)(1-q) & (1-p)q\\ p² & 2p(1-p) & (1-p)² \end{pmatrix} \]
Uma matriz de transição é válida se satisfaz as duas propriedades fundamentais:
\(P(x,y) \geq 0, \hspace{5mm} \forall x,y \in \mathcal{S}\).
\(\displaystyle \sum_{y \in \mathcal{S}}P(x,y) = 1, \hspace{5mm} x \in \mathcal{S}\).
\[ \mathcal{I} = \begin{pmatrix} + & + & +\\ + & + & +\\ + & + & + \end{pmatrix} \]
Uma cadeia de Markov de recorrência positiva e irredutível (como a nossa) possui uma única distribuição estacionária \(\pi\).
Além disso, satisfeitas as condições \(0 < p < 1\) e \(0 < q < 1\) a cadeia é aperiódica,e portanto ergódica. Logo, a distribuição estacionária existe e é igual a distribuição limite da cadeia.
Podemos obter essa distribuição em função de \(p\) e \(q\) ao resolver o sistema de equações resultantes de \(\pi = \pi P\).
\[ \begin{cases} (1 - q)^{2}\pi_{0} + p(1 - q)\pi_{1} + p^{2}\pi_{2} = \pi_{0} \\ 2q(1 - q)\pi_{0} + (pq + (1 - p)(1 - q))\pi_{1} + 2p(1 - p)\pi_{2} = \pi_{1} \\ q^{2}\pi_{0} + (1 - p)\pi_{1} + (1 - p)^{2}\pi_{2} = \pi_{2} \\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]
\[ \begin{cases} \pi_{0}[1−(1−q)^{2}] = \pi_{1}p(1−q) + \pi_{2}p^{2}\\ \pi_{2}[1−(1−p)^{2}] = \pi_{0}q^{2} + \pi_{1}(1−p)q \\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]
\[ \begin{cases} \pi_{0}[1−(1−q)^{2}] = \pi_{1}p(1−q) + \pi_{2}p^{2}\\ - \pi_{0}q^{2} = \pi_{1}(1−p)q - \pi_{2}[1−(1−p)^{2}]\\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]
\[ \begin{cases} A = 1−(1−q)^{2} = 2q−q^{2} \\ B = p(1−q) \\ C = p^{2} \\ D = 1−(1−p)^{2} = 2p−p^{2} \\ E = q(1-p) \\ F = q^{2} \end{cases} \]
\[ \begin{cases} B\pi_{1} + C\pi_{2} = A\pi_{0} \\ E\pi_{1} - D\pi_{2} = - F\pi_{0} \\ \pi_{0} + \pi_{1} + \pi_{2} = 1 \end{cases} \]
\[ \begin{cases} \pi_{0} = \frac{p^{2}}{(p+q)^{2}} \\[6pt] \pi_{1} = \frac{2q}{p} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{2pq}{(p+q)^{2}} \\[6pt] \pi_{2} = \frac{q^{2}}{p^{2}} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{q^{2}}{(p+q)^{2}} \end{cases} \]
Agora vamos realizar a implementação de uma simulação da cadeia desse sistema no R.
Partindo do estado 2, simularemos 1000 etapas da cadeia e iremos comparar as proporções dos tempos em que cadeia permance em cada estado com a distribuição estacionária.
Para a simulação que vem logo a seguir, os valores de \(p\) e \(q\) são escolhidos de forma arbitrária considerando os seguintes pontos:
Para a simulação que vem logo a seguir, os valores de \(p\) e \(q\) são escolhidos de forma arbitrária considerando os seguintes pontos:
library("tidyverse")
library("gt")
matrix_pot <- function(M, n){
# Verificacao de pressupostos.
# Matriz quadrada.
if(dim(M)[1] != dim(M)[2]){
print("Erro! A matriz nao quadrada.") # Avisa sobre o erro.
return(M * 0) # Retorna uma matriz nula para chamar atenção.
}
# Com entradas nao negativas.
if(any(M < 0)){
return("Erro! A matriz possui alguma entrada negativa.") # Avisa sobre o erro.
print(M * 0) # Retorna uma matriz nula para chamar atenção.
}
# A soma nas linhas deve ser igual a 1.
if(any(abs(apply(M, 1, sum) - 1) > 1e-09)){
print("Erro! Alguma linha na matriz nao soma 1.") # Avisa sobre o error.
return(M * 0) # Retorna uma matriz nula para chamar atenção.
}
x <- diag(dim(M)[1])
while(n){
x <- x %*% M
n <- n - 1
}
return(x)
}
simula_cadeia <- function(S, P, pi0, n){
# Simula uma cadeia de Markov.
#
# Parâmetros de entrada:
#
## S: espaco de estados finito de onde obteremos os estados indexados nas linhas e colunas da matriz de transição.
## P: matriz de transicao.
## pi0: vetor de probabilidades dos estados S.
## n: número de etapas a serem realizadas.
# Gerar X0 com distribuicao pi0 e armazenar no vetor x.
x <- numeric(n + 1)
x[1] <- sample(S,size = 1, prob = pi0)
# Gerar os estados das próximas etapas.
for(i in 1:n){
j <- match(x[i], S) # pega o índice do estado atual.
p <- P[j,] # pega a função de transição para o estado da etapa presente.
x[i+1] <- sample(S, size = 1, prob = p) # obter estado da proxima etapa.
}
return(x)
}
## Aplicação
# Já obtidos anteriormente.
p <- p # probabilidade de falha
q <- q # probabilidade de reparo
# Número de realizações
N <- 1000
P <- matrix(
c(
(1-q)**2, 2*q*(1-q), q**2,
p*(1-q), p*q + (1-p)*(1-q), (1-p)*q,
p**2, 2*p*(1-p), (1-p)**2
), byrow = TRUE, nrow = 3
)
# Fixar semente.
set.seed(123)
# Espaço de estados e distribuição inicial.
S <- c(0,1,2)
pi0 <- c(0,0,1)
# Simular a cadeia.
cadeia <- simula_cadeia(S, P, pi0, N)[-1]
cadeia_tbl <- tibble(etapa = 1:length(cadeia), estado = cadeia) |>
arrange(etapa)
tabela_resultados <- cadeia_tbl |>
group_by(estado) |>
count() |>
mutate(freq = n/N) |>
gt(groupname_col = FALSE) |>
tab_header(paste("Simulação: número de visitas a cada estado para", N, "etapas")) |>
cols_label(
estado = "Estado",
n = "Número de visitas",
freq = "Proporção de visitas"
)| Simulação: número de visitas a cada estado para 1000 etapas | ||
|---|---|---|
| Estado | Número de visitas | Proporção de visitas |
| 0 | 121 | 0.121 |
| 1 | 461 | 0.461 |
| 2 | 418 | 0.418 |
0 1 2
[1,] 0.1225 0.455 0.4225
[2,] 0.1225 0.455 0.4225
[3,] 0.1225 0.455 0.4225
\[ \begin{cases} \pi_{0} = \frac{p^{2}}{(p+q)^{2}} = \frac{0.35^{2}}{1^{2}} = 0.1225 \\ \pi_{1} = \frac{2q}{p} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{2pq}{(p+q)^{2}} = \frac{2 \cdot 0.35 \cdot 0.65}{1^{2}} = 0.4550\\ \pi_{2} = \frac{q^{2}}{p^{2}} \cdot \frac{p^{2}}{(p+q)^{2}} = \frac{q^{2}}{(p+q)^{2}} = \frac{0.65^{2}}{1^{2}} = 0.4225 \end{cases} \]
De acordo com os resultados da simulação anterior, a frequência de visitas aos estados em que o sistema está em operação (1 e 2) é bem maior que se comparada à frequência do estado em que o sistema está inoperante.
Isto se deve ao fato de que especificamos uma alta probabilidade do sistema se manter em funcionamento (65%) e uma baixa probabilidade do sistema permanecer inoperante (35%).
Os valores de \(p\) e \(q\) impactarão diretamente na qualidade e robustez do sistema como um todo.
Analisando-se uma máquina isolada (vide exemplo em aula de cadeia com dois estados), tem-se que a disponibilidade \(a\) de uma máquina no n-ésimo estado, isto é, \(P(X_{n} = 1)\) é dada por \(a = \frac{q}{p+q}\)
Ou seja, em regime estacionário, o número de máquinas funcionando segue uma distribuição \(Binomial(2, a)\), e portanto o número esperado de máquinas em funcionamento é \(2a\).
\[ P = \begin{pmatrix} 1 & 0 & 0\\ p & (1-p) & 0\\ p² & 2p(1-p) & (1-p)² \end{pmatrix} \]
\[ P = \begin{pmatrix} 1-p & 0 & p\\ 2p(1-p) & (1-p)^{2} & p^{2}\\ 0 & 0 & 1 \end{pmatrix} \]
\[ Q = \begin{pmatrix} 1-p & 0 \\ 2p(1-p) & (1-p)^{2} \end{pmatrix} \]
e
\[ R = \begin{pmatrix} p\\ p^{2} \end{pmatrix} \]
\[ (I-Q) = \begin{pmatrix} p & 0 \\ -2p(1-p) & p(2-p) \end{pmatrix} \]
\[ N = \begin{pmatrix} \frac{1}{p} & 0 \\ \frac{2(1-p)}{p(2-p)} & \frac{1}{p(2-p)} \end{pmatrix} \]
\[ R = \begin{pmatrix} \frac{1}{p}\\ \frac{3-2p}{p(2-p)} \end{pmatrix} \]
Observa-se que caso o sistema inicie com apenas uma máquina funcionando, o tempo médio de falha é geométrico \(\frac{1}{p}\);
Em ambos os casos, se \(p\) se aproxima de 1, o sistema tende a ficar inoperante em um passo;
Podemos analisar um grafico com os tempos médios de falha total para diferentes valores de \(p\), considerando-se o estado inicial 2.